home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zgeqpf.f < prev    next >
Text File  |  1996-07-19  |  7KB  |  224 lines
  1.       SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     March 31, 1993
  7. *
  8. *     .. Scalar Arguments ..
  9.       INTEGER            INFO, LDA, M, N
  10. *     ..
  11. *     .. Array Arguments ..
  12.       INTEGER            JPVT( * )
  13.       DOUBLE PRECISION   RWORK( * )
  14.       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  ZGEQPF computes a QR factorization with column pivoting of a
  21. *  complex M-by-N matrix A: A*P = Q*R.
  22. *
  23. *  Arguments
  24. *  =========
  25. *
  26. *  M       (input) INTEGER
  27. *          The number of rows of the matrix A. M >= 0.
  28. *
  29. *  N       (input) INTEGER
  30. *          The number of columns of the matrix A. N >= 0
  31. *
  32. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  33. *          On entry, the M-by-N matrix A.
  34. *          On exit, the upper triangle of the array contains the
  35. *          min(M,N)-by-N upper triangular matrix R; the elements
  36. *          below the diagonal, together with the array TAU,
  37. *          represent the orthogonal matrix Q as a product of
  38. *          min(m,n) elementary reflectors.
  39. *
  40. *  LDA     (input) INTEGER
  41. *          The leading dimension of the array A. LDA >= max(1,M).
  42. *
  43. *  JPVT    (input/output) INTEGER array, dimension (N)
  44. *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
  45. *          to the front of A*P (a leading column); if JPVT(i) = 0,
  46. *          the i-th column of A is a free column.
  47. *          On exit, if JPVT(i) = k, then the i-th column of A*P
  48. *          was the k-th column of A.
  49. *
  50. *  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
  51. *          The scalar factors of the elementary reflectors.
  52. *
  53. *  WORK    (workspace) COMPLEX*16 array, dimension (N)
  54. *
  55. *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
  56. *
  57. *  INFO    (output) INTEGER
  58. *          = 0:  successful exit
  59. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  60. *
  61. *  Further Details
  62. *  ===============
  63. *
  64. *  The matrix Q is represented as a product of elementary reflectors
  65. *
  66. *     Q = H(1) H(2) . . . H(n)
  67. *
  68. *  Each H(i) has the form
  69. *
  70. *     H = I - tau * v * v'
  71. *
  72. *  where tau is a complex scalar, and v is a complex vector with
  73. *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
  74. *
  75. *  The matrix P is represented in jpvt as follows: If
  76. *     jpvt(j) = i
  77. *  then the jth column of P is the ith canonical unit vector.
  78. *
  79. *  =====================================================================
  80. *
  81. *     .. Parameters ..
  82.       DOUBLE PRECISION   ZERO, ONE
  83.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  84. *     ..
  85. *     .. Local Scalars ..
  86.       INTEGER            I, ITEMP, J, MA, MN, PVT
  87.       DOUBLE PRECISION   TEMP, TEMP2
  88.       COMPLEX*16         AII
  89. *     ..
  90. *     .. External Subroutines ..
  91.       EXTERNAL           XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R
  92. *     ..
  93. *     .. Intrinsic Functions ..
  94.       INTRINSIC          ABS, DCMPLX, DCONJG, MAX, MIN, SQRT
  95. *     ..
  96. *     .. External Functions ..
  97.       INTEGER            IDAMAX
  98.       DOUBLE PRECISION   DZNRM2
  99.       EXTERNAL           IDAMAX, DZNRM2
  100. *     ..
  101. *     .. Executable Statements ..
  102. *
  103. *     Test the input arguments
  104. *
  105.       INFO = 0
  106.       IF( M.LT.0 ) THEN
  107.          INFO = -1
  108.       ELSE IF( N.LT.0 ) THEN
  109.          INFO = -2
  110.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  111.          INFO = -4
  112.       END IF
  113.       IF( INFO.NE.0 ) THEN
  114.          CALL XERBLA( 'ZGEQPF', -INFO )
  115.          RETURN
  116.       END IF
  117. *
  118.       MN = MIN( M, N )
  119. *
  120. *     Move initial columns up front
  121. *
  122.       ITEMP = 1
  123.       DO 10 I = 1, N
  124.          IF( JPVT( I ).NE.0 ) THEN
  125.             IF( I.NE.ITEMP ) THEN
  126.                CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
  127.                JPVT( I ) = JPVT( ITEMP )
  128.                JPVT( ITEMP ) = I
  129.             ELSE
  130.                JPVT( I ) = I
  131.             END IF
  132.             ITEMP = ITEMP + 1
  133.          ELSE
  134.             JPVT( I ) = I
  135.          END IF
  136.    10 CONTINUE
  137.       ITEMP = ITEMP - 1
  138. *
  139. *     Compute the QR factorization and update remaining columns
  140. *
  141.       IF( ITEMP.GT.0 ) THEN
  142.          MA = MIN( ITEMP, M )
  143.          CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
  144.          IF( MA.LT.N ) THEN
  145.             CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A,
  146.      $                   LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO )
  147.          END IF
  148.       END IF
  149. *
  150.       IF( ITEMP.LT.MN ) THEN
  151. *
  152. *        Initialize partial column norms. The first n elements of
  153. *        work store the exact column norms.
  154. *
  155.          DO 20 I = ITEMP + 1, N
  156.             RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
  157.             RWORK( N+I ) = RWORK( I )
  158.    20    CONTINUE
  159. *
  160. *        Compute factorization
  161. *
  162.          DO 40 I = ITEMP + 1, MN
  163. *
  164. *           Determine ith pivot column and swap if necessary
  165. *
  166.             PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 )
  167. *
  168.             IF( PVT.NE.I ) THEN
  169.                CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
  170.                ITEMP = JPVT( PVT )
  171.                JPVT( PVT ) = JPVT( I )
  172.                JPVT( I ) = ITEMP
  173.                RWORK( PVT ) = RWORK( I )
  174.                RWORK( N+PVT ) = RWORK( N+I )
  175.             END IF
  176. *
  177. *           Generate elementary reflector H(i)
  178. *
  179.             AII = A( I, I )
  180.             CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
  181.      $                   TAU( I ) )
  182.             A( I, I ) = AII
  183. *
  184.             IF( I.LT.N ) THEN
  185. *
  186. *              Apply H(i) to A(i:m,i+1:n) from the left
  187. *
  188.                AII = A( I, I )
  189.                A( I, I ) = DCMPLX( ONE )
  190.                CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1,
  191.      $                     DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK )
  192.                A( I, I ) = AII
  193.             END IF
  194. *
  195. *           Update partial column norms
  196. *
  197.             DO 30 J = I + 1, N
  198.                IF( RWORK( J ).NE.ZERO ) THEN
  199.                   TEMP = ONE - ( ABS( A( I, J ) ) / RWORK( J ) )**2
  200.                   TEMP = MAX( TEMP, ZERO )
  201.                   TEMP2 = ONE + 0.05D0*TEMP*
  202.      $                    ( RWORK( J ) / RWORK( N+J ) )**2
  203.                   IF( TEMP2.EQ.ONE ) THEN
  204.                      IF( M-I.GT.0 ) THEN
  205.                         RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 )
  206.                         RWORK( N+J ) = RWORK( J )
  207.                      ELSE
  208.                         RWORK( J ) = ZERO
  209.                         RWORK( N+J ) = ZERO
  210.                      END IF
  211.                   ELSE
  212.                      RWORK( J ) = RWORK( J )*SQRT( TEMP )
  213.                   END IF
  214.                END IF
  215.    30       CONTINUE
  216. *
  217.    40    CONTINUE
  218.       END IF
  219.       RETURN
  220. *
  221. *     End of ZGEQPF
  222. *
  223.       END
  224.